!-------------------------------------------------------------------------------------
subroutine diagnostics( dfu, Ez2D, Ex2D, Ey2D, f0vec, t)
!-------------------------------------------------------------------------------------

  use dimensions
  use mpi_wrap

  use fft_wrap
  use parameters, only : kz, pi, nz, runname, restart, SpeckleWidth
  use efield, only     : efield_get_epsilon, efield_get_E0z_center

  implicit none

  !-- global variables --

  complex, intent(in)    :: dfu(jv, jz, jy, jx, jDir)
  complex, intent(in)    :: Ez2D(jz,jy,jx), Ex2D(jz,jy,jx), Ey2D(jz,jy,jx)
  real, intent(in)       :: f0vec(jv)
  real, intent(in)       :: t


  !-- local variables --

  logical, save          :: first_call = .true.

  real       :: PhasePhi
  real, save :: PhasePhiOld
  real       :: KE, PotE, PotEx, PotEy, potEz, Phi0Energy
  real       :: nTrapped, fsq(1), fsq1(1)


  complex    :: EzFtInternal, EzFTtotal
            
  real       :: funiform(jv), funiform1(jv)

  complex    :: Phi1Dz(jz)
  real       :: E0zCenter(jz)

  integer    :: ix, iy, iz, iv, iDir  ! counters


  !-----------------------------------------------------------------------------------
       
       
      fsq1 = 0  
      do iv=1,jv
        funiform1(iv) = sum( dfu(iv,:,:,:,:) )
        fsq1 =  fsq1 + sum( ( f0vec(iv) + dfu(iv,:,:,:,:) )**2  )
      end do

      call mpi_reduce(funiform1, funiform, jv, &
           MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)

      call mpi_reduce(fsq1, fsq, 1, &
           MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)


      if (myid .ne. 0) return
     
      funiform = funiform /(jx*jy*jz)

      fsq = fsq*dz*dv 
      if (jx .ne. 1) fsq = fsq*dx
      if (jy .ne. 1) fsq = fsq*dy


      PotEx = 0.5*sum(abs(Ex2D(:,:,:))**2) / (jx*jy*jz)
      PotEy = 0.5*sum(abs(Ey2D(:,:,:))**2) / (jx*jy*jz)       
      PotEz = 0.5*sum(abs(Ez2D(:,:,:))**2) / (jx*jy*jz)

      potE = PotEx + PotEy + PotEz
        

      KE =  0.5 * sum( funiform(:) * vvec(:)**2 ) * dv

      call efield_get_E0z_center(E0zCenter)

      Phi1Dz(:) = Ez2D(:, jy/2+1, jx/2+1) - E0zCenter(:)
      call dfftw_execute_dft(plan_zforward, Phi1Dz, Phi1Dz)
      Phi1Dz=Phi1Dz/float(jz)
      EzFtInternal = Phi1Dz(1+nz)
        
      Phi1Dz(:) = Ez2D(:, jy/2+1, jx/2+1)
      call dfftw_execute_dft(plan_zforward, Phi1Dz, Phi1Dz)
      Phi1Dz=Phi1Dz/float(jz)
      EzFTtotal = Phi1Dz(1+nz)

      PhasePhi = -aimag(log(EzFtInternal)) ! really phase of E
        
      do while ( ( PhasePhi - PhasePhiOld ) <  - pi )
         PhasePhi = PhasePhi + 2*pi
      enddo
      do while ( ( PhasePhi - PhasePhiOld ) >  + pi )
         PhasePhi = PhasePhi - 2*pi
      enddo
        
      PhasePhiOld = PhasePhi   

     !------------------------------
       

      if (first_call .and. .not. restart) then

         open(3, file=trim(runname)//'.dat')

         write(3, "('% 1.t  2.PhiInternal  3.PhiToTal  4.-Phase(phi)  5.KineticEnergy')")
         write(3, "('% 6.Ex^2energy  7.Ey^2energy  8.Ez^2energy')")
         write(3,*)

      else

         open(3, file=trim(runname)//'.dat', position='append')

      end if

      if (first_call)  first_call = .false.

      PHi0Energy=0
      nTrapped=0.0

           write(3,501) t, 2*abs(EzFtInternal)/kz, 2*abs(EzFTtotal)/kz, &
           PhasePhi, KE, potEx, potEy, potEz

      close(3)

    501   format(  6400(e12.6,2x)  )


end subroutine diagnostics


!-------------------------------------------------------------------------------------
subroutine diagnostics_lin( dfu, Ex2D, Ez2D, uxMat, t)
!-------------------------------------------------------------------------------------

  use dimensions
  use mpi_wrap
  use fft_3d

  use parameters, only : runname, restart, ux0

  implicit none

  !-- global variables --

  complex, intent(in)    :: dfu(jv, jz, jy, jx, jDir)
  complex, intent(in)    :: Ex2D(jz,jy,jx)
  complex, intent(in)    :: Ez2D(jz,jy,jx)
  complex, intent(in)    :: uxMat(jz, jy, jx, jDir)
  real, intent(in)       :: t


  !-- local variables --

  complex    :: Ek(jz,jy,jx)
  real       :: Eksq(jz,jx)

  real       :: KE, PotEx, potEz
  real       :: funiform(jv), funiform1(jv)
  real       :: uxsq(jDir), uxsq1(jDir)

  integer    :: ix, iy, iz, iv, iDir

  logical, save          :: first_call = .true.
  integer, save          :: filenum = 0

  character*4            :: file_ext


  !-----------------------------------------------------------------------------------
       
       
      do iv=1,jv
         funiform1(iv) = sum( dfu(iv,:,:,:,:) )
      end do

      do iDir=1,jDir
         uxsq1(iDir) = sum ( ( uxMat(:,:,:,iDir) + ux0(iDir) )**2 )
      end do


      call mpi_reduce(funiform1, funiform, jv, &
           MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)

      call mpi_reduce(uxsq1, uxsq, jDir, &
           MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)

      if (myid .ne. 0) return
     
      funiform = funiform /(jx*jy*jz)
      uxsq     = uxsq /(jx*jy*jz)


      PotEx = 0.5*sum(abs(Ex2D(:,:,:))**2) / (jx*jy*jz)
      PotEz = 0.5*sum(abs(Ez2D(:,:,:))**2) / (jx*jy*jz)

      KE =  0.5 * sum( funiform(:) * vvec(:)**2 ) * dv

      !-- individual modes of Ex --

      Ek = Ex2D

      call fft_forward_z(Ek)
      call fft_forward_x(Ek)
      call fft_forward_y(Ek)

      Eksq = abs(Ek(:,1,:))
      Eksq = Eksq**2


      !------------------------------
       
      open(3, file=trim(runname)//'.dat', position='append')

      if (first_call .and. .not. restart) then
  
         write(3, "('% 1.t  2.KineticEnergy  3.Ex^2energy  4.Ez^2energy  5.Ux^2average')")
         write(3,*)

      end if

      write(3,501) t, KE, potEx, potEz, sum(uxsq)
 
      close(3)


      !-- individual modes

      filenum = filenum+1
 
!     open(5, file= trim(runname)//'.Ek', &
!          form='unformatted', access='direct',recl=jx*jz*8)
!     write(5, rec=filenum ) Eksq
!     close(5)
  

      !------------------------------

      first_call = .false.

    501   format(  6400(e12.6,2x)  )


end subroutine diagnostics_lin



!-------------------------------------------------------------------------------------
subroutine diagnostics_nlvp( dfu, Ez2D, t)
!-------------------------------------------------------------------------------------

  use dimensions
  use mpi_wrap


  use parameters, only : runname, restart, ux0

  implicit none

  !-- global variables --

  complex, intent(in)    :: dfu(jv, jz, jy, jx, jDir)
  complex, intent(in)    :: Ez2D(jz,jy,jx)
  real, intent(in)       :: t


  !-- local variables --

  real       :: KE, PotEx, potEz
  real       :: funiform(jv), funiform1(jv)

  integer    :: ix, iy, iz, iv, iDir

  logical, save          :: first_call = .true.
  integer, save          :: filenum = 0

  character*4            :: file_ext


  !-----------------------------------------------------------------------------------
       
       
      do iv=1,jv
         funiform1(iv) = sum( dfu(iv,:,:,:,:) )
      end do


      call mpi_reduce(funiform1, funiform, jv, &
           MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)


      if (myid .ne. 0) return
     
      funiform = funiform /(jx*jy*jz)

      PotEz = 0.5*sum(Ez2D(:,:,:)**2) / (jx*jy*jz)

      KE =  0.5 * sum( funiform(:) * vvec(:)**2 ) * dv


      !------------------------------
       
      open(3, file=trim(runname)//'.dat', position='append')

      if (first_call .and. .not. restart) then
  
         write(3, "('% 1.t  2.KineticEnergy  3.Ez^2energy')")
         write(3,*)

         write(*, "('% 1.t  2.KineticEnergy  3.Ez^2energy')")
         write(*,*)

      end if

      write(3,501) t, KE, potEz
      write(*,501) t, KE, potEz

      close(3)


      !------------------------------

      first_call = .false.

    501   format(  6400(e12.6,2x)  )


end subroutine diagnostics_nlvp


!--------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------


